import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
from scipy.interpolate import griddata
# 新增导入scipy.ndimage用于中值滤波
from scipy import ndimage


def read_gps_data(filename):
    """
    读取GPS数据文件
    假设数据格式为: longitude, latitude, elevation
    """

    try:
        # 尝试读取CSV格式数据
        data = pd.read_csv(filename)
        return data
    except:
        # 如果不是CSV格式，尝试按空格或制表符分割
        try:
            data = pd.read_csv(filename, delimiter='\s+')
            return data
        except:
            raise Exception("无法读取数据文件")


# 新增函数：平滑处理突变数值
def smooth_outliers(data, z_col='水底高程', threshold=2.0):
    """
    使用中值滤波平滑处理突变数值
    :param data: 输入数据
    :param z_col: 高程列名
    :param threshold: 异常值判断阈值（标准差倍数）
    :return: 处理后的数据
    """
    # 创建数据副本避免修改原始数据
    smoothed_data = data.copy()

    # 获取高程值
    z_values = smoothed_data[z_col].values

    # 计算均值和标准差
    mean_val = np.mean(z_values)
    std_val = np.std(z_values)

    # 识别异常值（超过均值±threshold*标准差的值）
    outlier_mask = np.abs(z_values - mean_val) > threshold * std_val

    # 对异常值应用中值滤波进行平滑处理
    # 创建一个临时数组用于滤波
    z_filtered = ndimage.median_filter(z_values, size=3)

    # 只将滤波后的值应用于异常点
    z_values[outlier_mask] = z_filtered[outlier_mask]

    # 更新数据
    smoothed_data[z_col] = z_values

    return smoothed_data


def calculate_elevation_distribution(data, z_col='水底高程'):
    """
    计算各0.1间隔值分布比例
    """
    # 获取高程数据并转换为0.1间隔的数值
    z_values = data[z_col].values
    z_01 = np.round(z_values / 0.1) * 0.1  # 转换为0.1间隔
    z_01 = np.round(z_01, 1)  # 确保精度为1位小数

    # 计算每个0.1间隔值的出现次数
    unique_values, counts = np.unique(z_01, return_counts=True)

    # 计算比例
    total_count = len(z_01)
    proportions = counts / total_count

    # 创建结果DataFrame
    distribution = pd.DataFrame({
        '高程值': unique_values,
        '出现次数': counts,
        '比例': proportions,
        '百分比': proportions * 100
    })

    return distribution


# 新增函数：计算两个数据集之间的体积差
def calculate_volume_difference(old_data,
                                new_data,
                                x_col='东坐标',
                                y_col='北坐标',
                                z_col='水底高程'):
    """
    计算两个数据集之间的体积差
    """
    # 计算重叠区域的坐标范围
    x_overlap_min = max(old_data[x_col].min(), new_data[x_col].min())
    x_overlap_max = min(old_data[x_col].max(), new_data[x_col].max())
    y_overlap_min = max(old_data[y_col].min(), new_data[y_col].min())
    y_overlap_max = min(old_data[y_col].max(), new_data[y_col].max())

    # 检查是否存在重叠区域
    if x_overlap_min >= x_overlap_max or y_overlap_min >= y_overlap_max:
        return 0, None, None, (x_overlap_min, x_overlap_max, y_overlap_min,
                               y_overlap_max)

    # 在重叠区域创建网格
    # xi = np.linspace(x_overlap_min, x_overlap_max, 100)
    # yi = np.linspace(y_overlap_min, y_overlap_max, 100)
    # Xi, Yi = np.meshgrid(xi, yi)

    # # 分别插值两个数据集
    # old_z = old_data[z_col].values
    # new_z = new_data[z_col].values
    # old_x, old_y = old_data[x_col].values, old_data[y_col].values
    # new_x, new_y = new_data[x_col].values, new_data[y_col].values

    # old_Zi = griddata((old_x, old_y), old_z, (Xi, Yi), method='cubic')
    # new_Zi = griddata((new_x, new_y), new_z, (Xi, Yi), method='cubic')

    # # 计算网格单元面积
    # dx = xi[1] - xi[0]
    # dy = yi[1] - yi[0]
    # cell_area = dx * dy

    # # 计算重叠区域内的体积差
    # mask = ~(np.isnan(old_Zi) | np.isnan(new_Zi))

    # # 修改: 计算带符号的体积差，以便区分增加和减少
    # volume_difference = np.sum((new_Zi[mask] - old_Zi[mask])) * cell_area

    # # 修改: 返回绝对值体积差和带符号的体积差
    # absolute_volume_difference = np.sum(
    #     np.abs(new_Zi[mask] - old_Zi[mask])) * cell_area

    return (x_overlap_min, x_overlap_max, y_overlap_min, y_overlap_max)


## 新增函数：生成两个数据集的对比图
def generate_comparison_plot(old_data,
                             new_data,
                             overlap_range,
                             x_col='东坐标',
                             y_col='北坐标',
                             z_col='水底高程'):
    """
    生成两个数据集在重叠区域的对比图
    """
    # 设置中文字体 - 改进版本
    plt.rcParams['font.sans-serif'] = [
        'SimHei', 'Microsoft YaHei', 'SimSun', 'FangSong', 'KaiTi'
    ]
    plt.rcParams['axes.unicode_minus'] = False

    x_min, x_max, y_min, y_max = overlap_range

    # 过滤数据，只保留在重叠区域内的数据点
    old_data_filtered = old_data[(old_data[x_col] >= x_min)
                                 & (old_data[x_col] <= x_max) &
                                 (old_data[y_col] >= y_min) &
                                 (old_data[y_col] <= y_max)]

    # 修复：使用new_data而不是old_data来过滤y_col
    new_data_filtered = new_data[(new_data[x_col] >= x_min)
                                 & (new_data[x_col] <= x_max) &
                                 (new_data[y_col] >= y_min) &
                                 (new_data[y_col] <= y_max)]

    # 在重叠区域创建网格
    xi = np.linspace(x_min, x_max, 100)
    yi = np.linspace(y_min, y_max, 100)
    Xi, Yi = np.meshgrid(xi, yi)

    # 分别插值两个数据集到重叠区域
    old_z = old_data_filtered[z_col].values
    new_z = new_data_filtered[z_col].values
    old_x, old_y = old_data_filtered[x_col].values, old_data_filtered[
        y_col].values
    new_x, new_y = new_data_filtered[x_col].values, new_data_filtered[
        y_col].values

    old_Zi = griddata((old_x, old_y), old_z, (Xi, Yi), method='cubic')
    new_Zi = griddata((new_x, new_y), new_z, (Xi, Yi), method='cubic')

    # 绘制对比图
    # 修改: 调整图形尺寸和布局以获得更好的显示效果
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))  # 修改为3个子图
    fig.subplots_adjust(wspace=0.3, right=0.95, left=0.05)

    # 确定统一的颜色映射范围，从0到30
    vmin, vmax = 0, 30

    # old数据图
    contour1 = axes[0].contour(Xi,
                               Yi,
                               old_Zi,
                               levels=15,
                               colors='black',
                               linewidths=0.5)
    contourf1 = axes[0].contourf(Xi,
                                 Yi,
                                 old_Zi,
                                 levels=15,
                                 cmap='coolwarm',
                                 vmin=vmin,
                                 vmax=vmax)
    axes[0].clabel(contour1, inline=True, fontsize=8)
    axes[0].set_xlabel('东坐标 (Easting)')
    axes[0].set_ylabel('北坐标 (Northing)')
    axes[0].set_title('old.csv 数据')
    # 修改: 调整颜色条大小和位置
    cbar1 = plt.colorbar(contourf1, ax=axes[0], shrink=0.8)
    cbar1.set_label('高程 (米)')

    # new数据图
    contour2 = axes[1].contour(Xi,
                               Yi,
                               new_Zi,
                               levels=15,
                               colors='black',
                               linewidths=0.5)
    contourf2 = axes[1].contourf(Xi,
                                 Yi,
                                 new_Zi,
                                 levels=15,
                                 cmap='coolwarm',
                                 vmin=vmin,
                                 vmax=vmax)
    axes[1].clabel(contour2, inline=True, fontsize=8)
    axes[1].set_xlabel('东坐标 (Easting)')
    axes[1].set_ylabel('北坐标 (Northing)')
    axes[1].set_title('new.csv 数据')
    # 修改: 调整颜色条大小和位置
    cbar2 = plt.colorbar(contourf2, ax=axes[1], shrink=0.8)
    cbar2.set_label('高程 (米)')

    # 差异图（只显示负值）
    diff_Zi = new_Zi - old_Zi
    # 将正值设为NaN，这样就只显示负值
    diff_Zi_positive_mask = diff_Zi > 0
    diff_Zi_modified = diff_Zi.copy()
    diff_Zi_modified[diff_Zi_positive_mask] = np.nan

    # 修改: 将等高线间隔从默认的1米改为0.2米
    diff_min = np.floor(np.nanmin(diff_Zi) * 5) / 5  # 确保包含所有数据的最小值
    diff_max = 0  # 最大值设为0，因为我们只关心负值
    diff_levels = np.arange(diff_min, diff_max + 0.2, 0.2)  # 0.2米间隔

    contour3 = axes[2].contour(Xi,
                               Yi,
                               diff_Zi_modified,
                               levels=diff_levels,
                               colors='black',
                               linewidths=0.5)
    # 修改: 自定义颜色映射，0点为浅红色，负值为深红色
    # 创建从深红色到浅红色的颜色映射
    from matplotlib.colors import LinearSegmentedColormap
    colors = ['darkred', 'lightcoral']
    cmap_black_red = LinearSegmentedColormap.from_list('red_gradient', colors)

    contourf3 = axes[2].contourf(Xi,
                                 Yi,
                                 diff_Zi_modified,
                                 levels=diff_levels,
                                 cmap=cmap_black_red)
    axes[2].clabel(contour3, inline=True, fontsize=8)
    axes[2].set_xlabel('东坐标 (Easting)')
    axes[2].set_ylabel('北坐标 (Northing)')
    axes[2].set_title('差异图 (new - old, 仅显示负值)')
    # 修改: 调整颜色条大小和位置
    cbar3 = plt.colorbar(contourf3, ax=axes[2], shrink=0.8)
    cbar3.set_label('高程差 (米)')

    plt.show()


# 新增函数：显示new.csv单独图像并支持手动绘制分割线
def show_new_data_with_interactive_line(new_data,
                                        x_col='东坐标',
                                        y_col='北坐标',
                                        z_col='水底高程'):
    """
    显示new.csv数据的单独图像，并支持手动绘制分割线以生成断面图
    """
    # 设置中文字体
    plt.rcParams['font.sans-serif'] = [
        'SimHei', 'Microsoft YaHei', 'SimSun', 'FangSong', 'KaiTi'
    ]
    plt.rcParams['axes.unicode_minus'] = False

    # 创建图形和轴
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
    fig.subplots_adjust(wspace=0.3)

    # 获取数据范围
    x_min, x_max = new_data[x_col].min(), new_data[x_col].max()
    y_min, y_max = new_data[y_col].min(), new_data[y_col].max()

    # 创建网格
    xi = np.linspace(x_min, x_max, 100)
    yi = np.linspace(y_min, y_max, 100)
    Xi, Yi = np.meshgrid(xi, yi)

    # 插值数据
    x_vals = new_data[x_col].values
    y_vals = new_data[y_col].values
    z_vals = new_data[z_col].values

    Zi = griddata((x_vals, y_vals), z_vals, (Xi, Yi), method='cubic')

    # 绘制new数据图
    contour = ax1.contour(Xi,
                          Yi,
                          Zi,
                          levels=15,
                          colors='black',
                          linewidths=0.5)
    contourf = ax1.contourf(Xi,
                            Yi,
                            Zi,
                            levels=15,
                            cmap='coolwarm',
                            vmin=0,
                            vmax=30)
    ax1.clabel(contour, inline=True, fontsize=8)
    ax1.set_xlabel('东坐标 (Easting)')
    ax1.set_ylabel('北坐标 (Northing)')
    ax1.set_title('new.csv 数据')
    cbar = plt.colorbar(contourf, ax=ax1, shrink=0.8)
    cbar.set_label('高程 (米)')

    # 添加交互功能说明
    ax1.text(0.05,
             0.95,
             '点击并拖动以绘制分割线\n双击结束绘制',
             transform=ax1.transAxes,
             fontsize=10,
             verticalalignment='top',
             bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

    # 存储用户绘制的线
    line_points = []
    line_obj = None

    def on_click(event):
        """处理鼠标点击事件"""
        if event.inaxes == ax1 and event.dblclick:
            # 双击结束绘制
            if len(line_points) >= 2:
                # 生成断面图
                generate_cross_section(new_data, line_points, ax2)
            return

        if event.inaxes == ax1 and event.button == 1:  # 左键单击
            # 添加点到线中
            line_points.append((event.xdata, event.ydata))

            # 更新显示的线
            update_line_display()

    def update_line_display():
        """更新显示的线"""
        global line_obj
        if line_obj:
            line_obj.remove()

        if len(line_points) >= 1:
            xs, ys = zip(*line_points)
            line_obj, = ax1.plot(xs,
                                 ys,
                                 'r-',
                                 linewidth=2,
                                 marker='o',
                                 markersize=4)

        plt.draw()

    def generate_cross_section(data, points, ax):
        """根据绘制的线生成断面图"""
        if len(points) < 2:
            return

        # 清除断面图轴
        ax.clear()

        # 计算沿线的点
        line_x, line_y = zip(*points)
        distances = [0]
        for i in range(1, len(line_x)):
            dx = line_x[i] - line_x[i - 1]
            dy = line_y[i] - line_y[i - 1]
            distances.append(distances[-1] + np.sqrt(dx * dx + dy * dy))

        # 在线上插值获取高程值
        elevations = []
        for x, y in points:
            # 找到最近的数据点
            dists = np.sqrt((data[x_col] - x)**2 + (data[y_col] - y)**2)
            closest_idx = np.argmin(dists)
            elevations.append(data.iloc[closest_idx][z_col])

        # 绘制断面图
        ax.plot(distances,
                elevations,
                'b-',
                linewidth=2,
                marker='o',
                markersize=4)
        ax.set_xlabel('沿线距离 (米)')
        ax.set_ylabel('高程 (米)')
        ax.set_title('断面图')
        ax.grid(True)

        plt.draw()

    # 连接事件处理函数
    fig.canvas.mpl_connect('button_press_event', on_click)

    plt.show()


# 修改: 新增函数计算所有差值为负的容积和
def calculate_negative_volume_difference(old_data,
                                         new_data,
                                         x_col='东坐标',
                                         y_col='北坐标',
                                         z_col='水底高程'):
    """
    计算两个数据集之间所有差值为负的容积和（即新数据低于旧数据的部分）
    """
    # 计算重叠区域的坐标范围
    x_overlap_min = max(old_data[x_col].min(), new_data[x_col].min())
    x_overlap_max = min(old_data[x_col].max(), new_data[x_col].max())
    y_overlap_min = max(old_data[y_col].min(), new_data[y_col].min())
    y_overlap_max = min(old_data[y_col].max(), new_data[y_col].max())

    # 检查是否存在重叠区域
    if x_overlap_min >= x_overlap_max or y_overlap_min >= y_overlap_max:
        return 0, None, None, (x_overlap_min, x_overlap_max, y_overlap_min,
                               y_overlap_max)

    # 在重叠区域创建网格
    xi = np.linspace(x_overlap_min, x_overlap_max, 100)
    yi = np.linspace(y_overlap_min, y_overlap_max, 100)
    Xi, Yi = np.meshgrid(xi, yi)

    # 分别插值两个数据集
    old_z = old_data[z_col].values
    new_z = new_data[z_col].values
    old_x, old_y = old_data[x_col].values, old_data[y_col].values
    new_x, new_y = new_data[x_col].values, new_data[y_col].values

    old_Zi = griddata((old_x, old_y), old_z, (Xi, Yi), method='cubic')
    new_Zi = griddata((new_x, new_y), new_z, (Xi, Yi), method='cubic')

    # 计算网格单元面积
    dx = xi[1] - xi[0]
    dy = yi[1] - yi[0]
    cell_area = dx * dy

    # 计算差值矩阵
    diff_Zi = new_Zi - old_Zi

    # 创建掩码，只考虑有效数据且差值为负的部分
    mask = ~(np.isnan(old_Zi) | np.isnan(new_Zi))
    negative_mask = (diff_Zi < 0) & mask

    # 计算所有差值为负的容积和
    negative_volume_sum = np.sum(np.abs(diff_Zi[negative_mask])) * cell_area

    return negative_volume_sum, old_Zi, new_Zi, (x_overlap_min, x_overlap_max,
                                                 y_overlap_min, y_overlap_max)


# 新增函数：查找最大正负差值
def find_max_positive_negative_differences(old_data,
                                           new_data,
                                           x_col='东坐标',
                                           y_col='北坐标',
                                           z_col='水底高程'):
    """
    查找两个数据集之间最大的正值差值和最大的负值差值
    """
    # 计算重叠区域的坐标范围
    x_overlap_min = max(old_data[x_col].min(), new_data[x_col].min())
    x_overlap_max = min(old_data[x_col].max(), new_data[x_col].max())
    y_overlap_min = max(old_data[y_col].min(), new_data[y_col].min())
    y_overlap_max = min(old_data[y_col].max(), new_data[y_col].max())

    # 检查是否存在重叠区域
    if x_overlap_min >= x_overlap_max or y_overlap_min >= y_overlap_max:
        return None, None, (x_overlap_min, x_overlap_max, y_overlap_min,
                            y_overlap_max)

    # 在重叠区域创建网格
    xi = np.linspace(x_overlap_min, x_overlap_max, 100)
    yi = np.linspace(y_overlap_min, y_overlap_max, 100)
    Xi, Yi = np.meshgrid(xi, yi)

    # 分别插值两个数据集
    old_z = old_data[z_col].values
    new_z = new_data[z_col].values
    old_x, old_y = old_data[x_col].values, old_data[y_col].values
    new_x, new_y = new_data[x_col].values, new_data[y_col].values

    old_Zi = griddata((old_x, old_y), old_z, (Xi, Yi), method='cubic')
    new_Zi = griddata((new_x, new_y), new_z, (Xi, Yi), method='cubic')

    # 计算差值矩阵
    diff_Zi = new_Zi - old_Zi

    # 添加: 剔除差值超过2米的数据
    diff_Zi[np.abs(diff_Zi) > 2.0] = np.nan

    # 创建掩码，只考虑有效数据
    mask = ~(np.isnan(old_Zi) | np.isnan(new_Zi) | np.isnan(diff_Zi))

    # 在有效数据中查找最大正值和最大负值
    if np.any(mask):
        max_positive_diff = np.nanmax(diff_Zi[mask])
        max_negative_diff = np.nanmin(diff_Zi[mask])
        return max_positive_diff, max_negative_diff, (x_overlap_min,
                                                      x_overlap_max,
                                                      y_overlap_min,
                                                      y_overlap_max)
    else:
        return None, None, (x_overlap_min, x_overlap_max, y_overlap_min,
                            y_overlap_max)


def main():
    """
    主函数
    """

    # 设置全局中文字体支持
    plt.rcParams['font.sans-serif'] = [
        'SimHei', 'Microsoft YaHei', 'SimSun', 'FangSong', 'KaiTi'
    ]
    plt.rcParams['axes.unicode_minus'] = False

    # 示例数据生成（如果没有实际数据文件）
    # 在实际使用中，请替换为您的数据文件路径
    try:
        # 读取GPS数据
        old_data = read_gps_data('old.csv')
        new_data = read_gps_data('new.csv')

        # 将列名映射为程序中使用的名称
        old_data.rename(columns={
            '东坐标': '东坐标',
            '北坐标': '北坐标',
            '水底高程': '水底高程'
        },
                        inplace=True)
        new_data.rename(columns={
            '东坐标': '东坐标',
            '北坐标': '北坐标',
            '水底高程': '水底高程'
        },
                        inplace=True)

        # 打印old和new的坐标范围
        print("old.csv 数据范围:")
        print(
            f"  东坐标范围: {old_data['东坐标'].min():.2f} - {old_data['东坐标'].max():.2f}"
        )
        print(
            f"  北坐标范围: {old_data['北坐标'].min():.2f} - {old_data['北坐标'].max():.2f}"
        )
        print(
            f"  水底高程范围: {old_data['水底高程'].min():.2f} - {old_data['水底高程'].max():.2f}"
        )

        print("\nnew.csv 数据范围:")
        print(
            f"  东坐标范围: {new_data['东坐标'].min():.2f} - {new_data['东坐标'].max():.2f}"
        )
        print(
            f"  北坐标范围: {new_data['北坐标'].min():.2f} - {new_data['北坐标'].max():.2f}"
        )
        print(
            f"  水底高程范围: {new_data['水底高程'].min():.2f} - {new_data['水底高程'].max():.2f}"
        )

        # 平滑处理突变数值
        old_data = smooth_outliers(old_data, '水底高程')
        new_data = smooth_outliers(new_data, '水底高程')
        old_data = old_data[old_data['水底高程'] >= 0]
        new_data = new_data[new_data['水底高程'] >= 0]
        # # 计算体积差
        # 修改: 接收新的返回值
        overlap_range = calculate_volume_difference(old_data, new_data)
        # 重新计算基于过滤后的数据
        dx = np.linspace(overlap_range[0],
                         overlap_range[1], 100)[1] - np.linspace(
                             overlap_range[0], overlap_range[1], 100)[0]
        dy = np.linspace(overlap_range[2],
                         overlap_range[3], 100)[1] - np.linspace(
                             overlap_range[2], overlap_range[3], 100)[0]
        cell_area = dx * dy

        # 新增: 计算所有差值为负的容积和
        negative_volume_sum, old_Zi_neg, new_Zi_neg, _ = calculate_negative_volume_difference(
            old_data, new_data)
        # 添加: 剔除差值超过2米的数据
        diff_Zi_neg = new_Zi_neg - old_Zi_neg
        diff_Zi_neg[np.abs(diff_Zi_neg) > 2.0] = np.nan
        mask_valid_neg = ~(np.isnan(old_Zi_neg) | np.isnan(new_Zi_neg)
                           | np.isnan(diff_Zi_neg))
        negative_mask = (diff_Zi_neg < 0) & mask_valid_neg
        negative_volume_sum = np.sum(np.abs(
            diff_Zi_neg[negative_mask])) * cell_area
        print(
            f"old.csv 和 new.csv 之间所有差值为负的容积和: {negative_volume_sum:.2f} 立方单位")

        # 新增: 查找最大正负差值
        max_pos_diff, max_neg_diff, _ = find_max_positive_negative_differences(
            old_data, new_data)
        if max_pos_diff is not None and max_neg_diff is not None:
            print(f"最大正值差值: {max_pos_diff:.2f} 米")
            print(f"最大负值差值: {max_neg_diff:.2f} 米")
        else:
            print("无法计算最大正负差值")

        # 生成对比图
        generate_comparison_plot(old_data, new_data, overlap_range)

        # # 显示new.csv的数据分布
        # distribution = calculate_elevation_distribution(new_data)
        # print("\nnew.csv 高程整数值分布比例:")
        # print(distribution.to_string(index=False))

    except Exception as e:
        print(f"文件未找到: {e}")

    # # 计算并显示各整数值分布比例
    # distribution = calculate_elevation_distribution(new_data)
    # print("\n高程整数值分布比例:")
    # print(distribution.to_string(index=False))


if __name__ == "__main__":
    main()
